Mesoscopic lattice Boltzmann modeling of soft-glassy systems: theory and simulations 
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A multi-component lattice Boltzmann model recently introduced (R. Benzi et at, Phys. Rev. Lett. 
102, 026002 (2009)) to describe some dynamical behaviors of soft-flowing materials is theoretically 
analyzed. Equilibrium and transport properties are derived within the framework of a continuum 
free-energy formulation, and checked against numerical simulations. Due to the competition between 
' short-range inter-species repulsion and mid-range intra-species attraction, the model is shown to 

\ give rise to a very rich configurational dynamics of the density field, exhibiting numerous features 

, of soft-flowing materials, such as long-time relaxation due to caging effects, enhanced viscosity 

04 ' and structural arrest, ageing under moderate shear and shear-thinning flow above a critical shear 

,__( ' threshold. 

^3 : 

(3 ■ I. INTRODUCTION 

m '. 

' The study of the rheology of flowing soft systems, such as emulsions, foams, gels, slurries, colloidal glasses and 
f-j— ^ • related complex fluids, is gaining an increasing role in modern science and engineering [l], [3, H, |3j S B S 111 • From the 
Q theoretical standpoint, much of the fascination of these systems stems from the fact that they do not fall within any 
c/2 of three basic states of matter, gas-liquid-solid, but live rather on a moving border between them. Foams are typically 
a mixture of gas and liquids, whose properties can change dramatically with the changing proportion of the two; wet- 
foams can flow almost like a liquid, whereas dry-foams may conform to regular patterns, exhibiting solid-like behavior 
0. Emulsions can be paralleled to bi-liquid foams, with the minority species dispersed in the dominant (continuous) 
one. The behavior and, to same extent, the very existence itself of both foams and emulsions are vitally dependent on 
(— I surface tension, namely the interactions that control the physics at the interface between different phases/ components. 
Q ' Indeed, the presence of surfactants, i.e. a third constituent with the capability of lowering surface tension, has a 
O ; profound impact on the behavior of foams and emulsions; by lowering the surface tension, surfactants can greatly 
facilitate mixing, a much sought-for property in countless practical endeavors, from oil-recovery, to chemical and 
biological applications. Another basic property of foams and emulsions is metastability /disorder. Indeed, in most 
^ ' instances, these materials consist of a disordered collection of droplets/bubbles with a broad distribution of sizes, 
randomly mixed and arranged, which do not correspond to the (global) minimum of any thermodynamic function. 
This is even truer in the case of complex flowing systems, which live consistently out of (thermodynamic) equilibrium. 
As a result, they exhibit a number of distinctive features, such as long-time relaxation, anomalous viscosity, aging 
. . behav ior, whose quantitative description is calling for profound extensions of non-equilibrium statistical mechanics 
' P, [13, [Hi [3 d, [13, EE [H, [13, [Hi ■ The study of these phenomena sets a pressing challenge for computer simulation 
. as well, since characteristic time-lengths of disordered fluids can escalate tens of decades over the molecular time scales 
' [l^ [13, [HI ■ In addition, tracking the time evolution of complex interfaces represents a serious hurdle for traditional 
' discretization techniques. These split into two broad categories: Eulerian and Lagrangian. In Eulerian methods, 
the physical observables are attached to a fixed grid and monitored as they change in time at each grid location. 
Lagrangian methods, on the contrary, "go with the flow", i.e. the degrees of freedom are attached to the moving 
fields, and most notably to the critical regions of the flow where the most abrupt changes take place (interfaces). 
cd As usual, both methods have their merits and pitfalls. Lagrangian methods do not waste degrees of freedom on 
uninteresting regions of the flow; however, since the grid adapts to the changing fields, when these changes are too 
abrupt the numerics is forced to ad-hoc readjustments (grid-rezoning) which may eventually fail and lead to collapse 
of the numerics [l^l . Eulerian methods are free from these problems, because the interface is not tracked, but just 
tagged as the region where strong gradients are detected. The downside is that very high resolution is needed around 
the interface, for otherwise excessive smoothing (numerical diffusion) results (diffuse-interface) [2^ . A special variant 
of Eulerian methods does not attempt to resolve the interface, which is treated as a zero-thickness mathematical 
interface, across which jumps of the observables are specified. A proper handling of the discontinuities, and the 
avoidance of spurious oscillations, is however a non- trivial task (24j . 

On a more microscopic scale, often too small for hydrodynamic purposes, to date the most credited techniques for 
complex fiowing materials are Molecular Dynamics and Monte Carlo simulations [l^ [2^, [2l| . Molecular Dynamics 
in principle provides a fully ab-initio description of the system, but it is limited to space-time scales significantly 
shorter than experimental ones. Monte Carlo methods are somehow less affected by this limitation, since they can be 
designed in compliance with accelerated-dynamic sampling rules. However, these rules meet with some difficulties in 
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FIG. 1: Sketch of the interaction forces between two species, say A and B. The two components A and B interact via a repulsive 
pseudo-potential, which supports a surface tension uab- Moreover, each component experiences an attractive interaction in 
the first Brillouin zone and a repulsive one acting on both Brillouin zones (see also figure [2] for the technical details). Each of 
these interactions can be tuned through a separate coupling constant. 

accounting for hydrodynamic interactions [1^. As a result, neither MD nor MC can easily take into account the non- 
equilibrium dynamics of complex flowing materials, such as micro-emulsions, on space-time scales of hydrodynamic 
interest. Besides these general techniques, a number of specialized methods are also available, such as dissipative 
particle dynamics [1^ and others. More in details, a special kind of Molecular Dynamics (MD) for dry granular, 
Stokesian Dynamics (SD) for viscous suspensions and Bubble- Model for foams [13, [H, [2^ have been developed with 
several adjustable particle interactions in order to have a good agreement with experiments [10, [3l|, [12] . By contrast, 
the combination of both particle deformation and viscous flow has not been fully described yet, although it is central 
in such materials as foams and emulsions. For macroscopic complex flows, particularly interesting is the approach by 
Doi et al. [33] , which construct a set of evolution equations for the volurne fraction of the oriented interface elements 
within a complex flows, and more recently the soft dynamics approach [34j . 

In the last decade, a new class of mesoscopic methods, based on minimal lattice formulations of Boltzmann's 
kinetic equation, have captured signiflcant interest as an efficient alternative to continuum methods based on the 
discretization of the Navier-Stokes equations for non-ideal fluids [H, [13, [H, [13, [HI . A very popular mesoscopic 
technique is the pseudo-potential-Latticc-Boltzmann (LB) method, developed over a decade ago by Shan & Chen 
[43 , 143I ]. In the SC method, potential energy interactions are represented through a density-dependent mean-field 
pseudo-potential, ^'(/o), and phase separation is achieved by imposing a short-range attraction between the light and 
dense phases. In this work, we discuss extensions of two-species, mesoscopic lattice Boltzmann model which prove 
capable of reproducing some features of fiowing soft- materials, such as structural arrest, anomalous viscosity, cage- 
effects and ageing under shear [3]. The key feature of the model is the capability to investigate the rheology of 
these systems on space-time scales of hydrodynamic interest at an affordable computational cost. Among others, this 
model shows the first evidence of mesoscopic cage formation and rupture within a hydrodynamic lattice Boltzmann 
description. 

The present work is organized in two major parts: Theory and Numerical Results. In section II, we provide 
the basic elements of the multicomponent lattice kinetic model with multi-range non-ideal interactions, short-range 
attraction and mid-range repulsion. In section III, we derive the macroscopic equations associated with the large-scale 
hydrodynamic limit of the kinetic model. In section IV and V, we present an explicit calculation of the equilibrium 
(equation of state) and transport (surface-tension) properties, both for the case of intra-species repulsion alone, as well 
as its combination with intra-species attraction. In the process, we detail how the combination of this short/mid-range 
attractive/repulsive interactions allows to bring the surface tension down to vanishingly small values, a property which 
is key to the complex and heterogeneous dynamics displayed by the model, and notably by the density field. The 
numerical part follows in section VI. In section VII we discuss the morphological features of the density configurations, 
and demonstrate the existence of long-lived metastable states resulting from the interplay/competition between short- 
range attraction and mid-range repulsion. In section VIII, we investigate the dynamic response of system under an 
external shear drive, and provide several evidences of complex behaviors, such as cage formation and rupture under 
shear, ageing and its disappearance above a critical shear threshold, long-term non-Newtonian shear-strain correlations 
and Barkhausen intermittcncy, namely a power-law distribution of the waiting times between sliding events events. In 
section IX we discuss the issues of sensitivity to initial conditions and finite-size effects. In section X, we conclude with 
an outlook and future perspectives for the application of the present model, and generalizations thereof, to a broad 
class of complex soft-flowing systems, such as foams, emulsions and similar. Finally, in the Appendix we provide the 
conversion rules from/to lattice to physical units. 
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II. THE MULTI-COMPONENT KINETIC MODEL 



Kinetic theory and its discrete (lattice Boltzmann) coun terp arts for multicomponent fluids and gas mixtures have 
received much attention in the literature [13, 111, 111, 111, [ij 111, [11 [11 [Sll [11 Many of the kinetic models 
for mixtures are based on the linearized Boltzmann equations, especially the single-relaxation-time model due to 
Bhatnagar, Gross, and Krook -the celebrated BGK model [s^. Here we shall consider the multicomponent model 
introduced by Shan & Chen [H [1^: after a brief summary the main properties of the model we will proceed 
to analyze the equilibrium states relevant on the hydrodynamic scales. We start from a kinetic lattice Boltzmann 
equation [3l,[43,[4l[ for a multicomponent fluid with Ng species [13, [ll] whose evolution equations over a characteristic 
time lapse At read as follows 

where fisi^, t) is the probability density function of finding a particle of species s = l-.-N^ at site r and time t, moving 
along the i-th lattice direction defined by the discrete speeds c; with i = 0...b. For simplicity, the characteristic time 
lapse At is assumed to be equal to unity in the following. The left hand-side of (|T]) stands for molecular free-streaming, 
whereas the right-hand side represents the time relaxation (due to collisions) towards local Maxwellian equilibrium 

fi's'^\psiu) on a time scale Ts [H, [H, [iO, [iH . The local Maxwellian is truncated at second order, an approximation 
that is sufficient to recover correct hydrodynamic balance in the isothermal regime 

fis {Ps,U) = Wl ^'ps 1 + 2 + TTir^ ^aUb 

with c| the square of the sound speed velocity in the model and Sab the Kronecker delta with a, b indicating the 

Cartesian components (repeated indices are summed upon). The w['"^^^s are equilibrium weights used to enforce 
isotropy of the hydrodynamic equations [ssl . [4oL [4lj | . To be noted that the equilibrium for the s species is a function 
of the local species density 



,(f,t) -^^(f,i) 



and the common velocity defined as 



u(f, t) 



This common velocity receives a shift from the force Fg acting on the s species [43 . [54| . This force may be an external 
one or it could also be due to intermolecular (pseudo)-potential interactions. The pseudo-potential force within each 
species consists of an attractive (a) component , acting only on the first Brillouin region (belt, for simplicity), and a 
repulsive (r) one acting on both belts, whereas the force between species (X) is short-ranged and repulsive: 



F,(r, t) = F:{r, t) + Fl{r, t) + F^ {f, t) 



where 



i£beltl 

Fl{r,t) - -Gl^s{r,t) ^ pS! ~ G^^ ,{r,t) ^ pS's{n,t)c, (2) 



iGbeltl iebelt2 



(Po '') s'^s iebeltl 

In the above, the groups 'belt 1' and 'belt 2' refer to the first and second Brillouin zones in the lattice and Ci, pi, Wi 
are the corresponding discrete speeds and associated weights (see figure [2] and tabic [l| . Apart from a normalization 
factor, these correspond to the values given in [1^, [13|. Also, Ggs' ~ Gs's, s' ^ s, is the cross-coupling between 
species, po a reference density to be defined shortly and, finally, fi — f + Ci are the displacements along the Ci 
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velocity vector. These interactions are sketched in Figure [T] for the case of a two component fluid (say species A 
and B). Note that positive (negative) G code for repulsion (attraction) respectively. This model is reminiscent of the 
potentials used to investigate arrested phase-separation and structural arrest in charged-colloidal s^^stems, and also 
bears similarities to the NNN (next-to-nearest-neighbor) frustrated lattice spin models [H, [H, [13, [6l[ . As compared 
with lattice spin models, in our case a high lattice connectivity is required to ensure compliance with macroscopic 
non-ideal hydrodynamics, particularly the isotropy of potential energy interactions, which lies at the heart of the 
complex rheology to be discussed in this work. To this purpose, the first belt is discretized with 9 speeds, while the 
second with 16, for a total of 6 = 25 connections (including rest-particles, for normalization purposes). The weights 
are chosen in such a way as to fulfill the following normalization constraints: 

Wo + ^ = po + ^ K + ^ Pi = 1 (3) 

iebeltl iGbeltl i£belt2 

= H p^'^^^^ + P'^^^^ = ^^s- (4) 

iebeltl iebeltl i£belt2 

with c| = 1/3 the lattice sound speed and wq and po the weights associated to the velocity at rest. All the weights 
take the values illustrated in Table [H The set of discrete speeds and corresponding weights are such as to recover 
4th order isotropy for the interactions running on the first belt and 8th order isotropy for those extending over the 
second one. This choice is naturally patterned after reference [l^ll^], although differen t op tions might be available. 
The pseudo-potential ^s{ps) is taken in the form originally suggested by Shan & Chen [42ll43| 

^[ps)^p'i\l~e-P^/''o'), (5) 

where marks the density value above which non ideal-effects come into play for species s. For the sake of simplicity, 
in the sequel we shall take a common value for all species, pj, = Po- 



Forcing Weigths (for FJ) 



Pi = 247/420 


i = 


P. = 4/63 


i= 1,4 


p» = 4/135 


i = 5,8 


p» = 1/180 


i = 9, 12 


p» = 2/945 


i = 13, 20 


Pi = 1/15120 


i = 21,24 



Forcing Weights (for and F^^) 


w, = 4/9 


i = 


10, = 1/9 


j = 1,4 


Wi = 1/36 


j = 5,8 



TABLE I: Links and weights of the two belts, 25-speeds lattice [stI. [6^ for all interactions sketched in equations ([2}. The first 
belt lattice velocities are indicated with i = 1...8 while the second belt ones with i = 9. ..24 (see also figure[2]for a sketch), pi or 
Wi is indicating the weight associated with the i-th velocity in the various interactions. The weights associated to the velocity 
at rest, wq and po, are chosen to enforce a unitary normalization Q. 



III. MACROSCOPIC EQUATIONS 

The set of macroscopic equations associated with our kinetic model consists of the continuity equations, one for 
each component separately, plus an equation of motion for total fluid momentum. Under the assumption of the same 
characteristic time scale for all the components Ts = r, [s^ these equations read as follows: 



dtPs + daiPsUa) = daJsa 



(6) 



23 



18 12 



19 



24 



FIG. 2: The discrete 25-speed lattice. Both behs are iUustrated with the corresponding velocities. 



dt{pUa) + dbipUaUb) = -db{clp + (Tab) + ^sa = -db{Pab - f^ab ) C^) 



where p = Ps is the total density, u — J2s Ps^s/p is the baricentric (total) fluid velocity, Fga the a-th component 

of the force acting on spec 
current in (|18p is given by 



of the force acting on specie s and cr^"^'^"' the dissipative component of the momentum-flux tensor. The diffusive 



Jsa=cl[T--] {daPs~-dap] ~T{F,a~-y^F,,a] . (8) 



Central to this analysis is the momentum-flux tensor, defined as the sum of a kinetic component plus an interaction 
term: 

Pab = + P'aT (9) 

where 

j^ar = E/-^«c., (10) 

is 

plus the interaction component, -P^^*, defined by the condition: 

Taylor expansion of the forcing terms will allow for a direct computation of P^"* and the diffusion currents [sj, [s^ . 
It has to be noted that relation pT|) can also be directly satisfied on the lattice using the idea developed in a recent 
paper by Shan [63| . thus leading to more refined computational results for the momentum equation. 



A. Two component fluid 

The picture simplifies significantly for the case of a two-component fluid (say A and B). When the distribution 
functions fiA, fiB are close to the equilibrium, the kinetic part of the pressure tensor takes the following form: 

Pai," = iPA+PB)45ab + KV (12) 
^(t) ^ ^4 PAPB f ^aPA _ 9qPb ^ ^ dbPA _ ^AB ^ ^^^^ 



J \ PA PB J \ PA PB 
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where we recognize cin ide^l ps-rt (^pA ~l~ PB^^g^ab 

plus some extra r dependent terms. The origin of these terms will 
be elucidated in the section devoted to transport properties when we will detail the calculations of the surface tension 
coefficients across curved interfaces. Upon Taylor expanding [s^ up to the fourth order the forcing terms in the 
momentum equation, the interaction terms of the pressure tensor can be recast into the following form 

PIT = (4%i*i + 4%i*l + clgABPAPB + 4n) s,, - r,, (m) 

n= E Gs2(\\V^s? + \^s^^)\ +^^(fiA^PB+ Pb^Pa + Vpa-^Pb) (15) 



s=A,B 



Tab = ^4 {G2Ada^Adb^A + G^Bda^ Bdb^ B + gAB{daPAdbPB + daPBdbPA)) ■ (16) 

In the above, we have set 

Qab = Gab/Po 

and introduced the effective couplings 

Gsi^g: + g: Gs2^g: + ^g: s^a,b. (i?) 

From these expressions, we note that the presence of the second-neighbor repulsive layer allows a separate control 
of the equilibrium (equation of state, i.e. terms proportional to c|) and transport properties (surface tension, i.e. 
terms proportional to Cg). For the diffusive current, we can Taylor expand the forcing terms up to the second order 
to obtain 

Jsa^^Dss'{pA,PB)daPs' S,s'^A,B (18) 

s' 

where the (non-Hnear) diffusion coefRcients are given by: 

Daa = c| (^y " ^) + J(G^i/^s*A*A - OabPePa)^ (19) 

Dbb = 4 (^y " + ^(Gbi/9a*b*b - 9abPaPb)^ (20) 

These are nothing but equations (26)- (29), already discussed in a earlier paper by Shan & Doolen js^l- The above 
expressions indicate that the intra-species mass flow consists of an internal component, proportional to the density 
of the other species, and a force-induced component, proportional to the intermolecular couplings [13, H^. Note that 
the latter does not vanish even in the limit of zero inter-species interactions, qab ^ 0. The following reciprocity 
relations: 

Dab = ~Dbb, Dba = -Daa- (21) 
secure conservation of the total density. The continuum-time limit t 3> ^ is thus characterized by 

Daa cIt f ^ + -{GaiPb-^ a'^'a " 9abPbPa)] (22) 
\ P P J 



Dbb c^T [ — + -{GbiPa'^ b'^'b - 9abPaPb)\ 
\ P P J 



(23) 



with the relaxation properties factorizing outside. It is therefore natural and convenient to introduce a r-dependent 
parameter 



(r) = ^ (24) 
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measuring the importance of discrete-time effects in the macroscopic equations. Clearly, in the continuum time limit 
6(t ) — > 1, while for r ^ 1/2 we have 9{t) — (in terms of lattice Boltzmann fluids this is a dissipation free limit 
plj). The diffusion coefficients can thus be recast into the following form 

Daa = rcl (e{T)^ + -{GaiPb-^a^'a ^ 9abPbPa)] (25) 



Dbb =tcI[0{t)^ + - (Gbipa*b*b - QabPaPb) ] ■ (26) 
P P 



IV. EQUILIBRIUM PROPERTIES 

In this section we study the main equilibrium properties of the model previously introduced whenever stable inter- 
faces between the two fluids set in. To this purpose, we will focus on a one-dimensional problem, where inhomogeneities 
in the density profiles develop only across a single coordinate, say x. It proves expedient to start with the case of two 
components with mutual density repulsion (i.e. equation ^ with Fl = 0, = ), where an exact matching with 
a free-energy functional can be achieved in the continuum limit, i.e. when the discrete lattice effects are negligible. 
This allows us to envisage efhcient strategies to describe the bulk equilibrium properties in special situations where 
all pseudo-potentials interactions are included (i.e. equation ([2]) with all the interactions on). 



A. Multicomponent Model with pure Density Repulsion 

At equilibrium, the relevant properties of the interfaces emerging from the separation of the fluids can be obtained 
by imposing a constant diffusion current and a constant pressure all across the interface (zero net flow can safely be 
assumed). This yields: 

rc| \ —0{t) - -{gABPBPA)] dxPA - tc| ( —0{t) - -{gABPAPB)) d^PB = Jo (27) 



P P J \ P P 

P'xx = c%pA + cIpb + clgABPAPB + Cs^^ [paOxxPb + PbQxxPA - dxPAdxPB) + K^jJ = Po (28) 

where Pq is the constant value of the pressure across the interface and Jq is the constant diffusion current predicted by 
the single component continuity equation. For simplicity we have not expanded the extra r dependent terms (KxV) of 
the kinetic pressure tensor ([TH)) . Since Jq = in the bulk phases {dxPA,B = 0), one concludes that Jg = everywhere. 
Next, we observe that the equation P7|) can be recast in the form of a differential equation relating the values of the 
two densities at each spatial location: 

dp A pAPg - PAPB 



dPB ORO^n^ 



In the above, we have defined 



PBPg - PAPB 



P'P - ^ (29) 

dAB 



as a characteristic density depending both on the relaxation properties in Bij) and on the intermolecular coupling 
gAB-i above which inter-species repulsion becomes dominant. At the spatial location where pA = PB^ we also have 
dxPA = —dxPB because of the symmetry of the system upon the interchange pA ^ Pb- Equation ([27| also shows that, 

(r) 

at this location, pA = PB = Pg ■ By integrating the previous differential equation backward and forward in density 
space, starting from the point where = —1, it is possible to construct the manifold of density pairs (pA, Pb) 
obeying the condition of zero mass flow. For the specific case in point, these equations can be solved exactly, leading 
to the following relation 



PA_ 

PB 



exp (^{pA - Pb)/p''J'^^ 
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Obviously, this relation is fulfilled by the trivial solution pa = Pb\ owing to the non-linearity of the above equations, 
non trivial solutions are expected beyond a critical value of Qab- These identify with the bulk densities once separation 
between the fluids has occurred. 

Since we have neglected higher order terms in the Taylor-expansion yielding the diffusive current, this relationship 
is not expected to hold uniformly across the interface. However, it can be be regarded as an excellent approximation 
to compute the bulk densities after separation of the two fluids. To this end, we note that, out of the full set of 
pair densities, {pA, Pb) belonging to the density manifold, only one is compatible with the condition of equilibrium. 
Mechanical equilibrium, as obtained by imposing a constant pressure tensor across the interface, equation (j28p . cannot 
serve as a selection criteria, because of the invariance under the interchange pA Pb- The two values of the bulk 
densities can however be fixed by imposing the total average density in the numerical simulations {pA + Pb) = {p)- 
This provides a system of two equations determining the two bulk densities: 

\^-^-p{{pa-Pb)Ip'^') (3P^ 

\PA + PB = {p)- 

Once the bulk densities have been fixed, the momentum equation (j28p . consistently with the higher order in the 
Taylor expansion for the density equation ([27|) . would allow to reconstruct the profiles across the interface. Such 
technical construction will make the object of a forthcoming paper. 



B. Free-energy procedure 



In order to better elucidate the mechanism fixing the bulk densities in the phase separation process, we can also 
resort to a direct exact link with a free energy functional in the continuum limit, where all discrete lattice effects 
disappear. We begin by considering a free-energy density in the form 



^PA,Pb) = clpA^OgPA + clpB^OgpB + clgABPAPB " Cg^^V pA ' Vpfi. 



(31) 



This consists of the sum of two ideal free-energy densities {cgPA,B^og pa.b) plus an interaction term. It has to be 
stressed that the terms proportional to CgQAB in front of the interacting terms should by no means be related to the 
fluid temperature, as they simply disappear upon a suitable choice of the lattice forcing weights [6^. On the other 
hand, the term proportional to c| in front of the ideal parts ('^ PA,B^ogpA,B) plays the role of a global reference 
temperature. This is of no relevance for the present athermal case, but may become important for generalizations 
involving temperature fluctuations [6^. where internal energies need to be introduced. As to the free-energy in (j3ip . 
it is readily checked that the bulk contribution 



fb{PA,PB) 

correctly reproduces the bulk pressure: 



C%pA log PA + C^sPB log /OS + clgABPAPB 



U / \ 9fb Ofb 2 t , \ , 2 

^b(PA,PB) = PA^ h PBt: /fc = Cg{pA + Pb)+ CgQABPAPB 

OpA OpB 

that is the generalization of the standard Legendre's relation Pb{p) = P^^f^ — f{p) connecting the free-energy to 
the bulk pressure of a single-component fluid. In order to preserve both densities separately, we next introduce two 
Lagrange multipliers, say Xa and As, thus leading to the following constrained free-energy density: 



^PA,Pb) = fb{pA,PB) - Cg^^VpA ' Vps - XaPA - ^bPB- 



(32) 



Variations of this constrained free-energy with respect to pA and ps delivers the following two Euler-Lagrange equa- 
tions: 







( "\ 


j dpA 






- da 


( 1 


y dpB 





Based on ([31]) . these yield: 



dp A ^^S 2 OxxPB 
^ + C%^dxxPA 



= 

= 0. 



As. 



(33) 



(34) 
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Upon multiplying the first equation by and the second by we can then integrate between the bulk region 
(x = 0) and a generic interface location (x). In this way, we obtain 

- ^sPa\1 - clgABPAPB\l + C%gAB Jo PBdyPAdy + 4^ Io{dyPA){dyyPB)dy - ci;^pAd^j:PB = ^g^^ 

- cIpb\1 - clgABPAPB\l + c^s9AB Jq pAdyPBdy + {dyPB){dyyPA)dy - c%^pBdxxPA = 

where ... |q represents the variation between (bulk) and x (interface location) of the desired observable. The above 
equations represent the conserved currents associated with the two Lagrange multipliers and they can be linked 
directly into the constant pressure tensor and diffusion current at equilibrium. In fact, by summing both equations 
in (|35p we obtain 

c%{pA + Pb) + clgABPAPB + '^^^^^ (pbOxxPa + PaOxxPb - dxPAdxPs) = const. (36) 

that is reminiscent of ([^5]) upon neglecting KxV . Similarly, upon applying the x derivative to both equations in ([55]) 
and then multiplying the first equation by pB and the second by pA we can finally subtract the two contributions to 
get 

{PbOxPA - PaQxPb) " 9ABPBPA{dxPA " QxPb) + 0(9^) = 

that is delivering the condition of a zero diffusion current, as given in (|27p with Jq = 0, in the limit r ^ i, i.e. 
9{t) 1. Thus, in the free-energy formalism, both conservations descend from the same single scalar. The free- 
energy formalism permits to recast the continuity equations in terms of the gradients of the chemical potentials 
Ps = More specifically: 

dtPA + daifAUa) = da{M {p A, PB)da{p.A - Pb)) (37) 



dtPB + daipBUa) = da{M {pA, PB)da{pB ~ Pa)) (38) 

where the mobility M{pA, Pb) is given by AI{pA, Pb) = r P^P^ , The above form of the continuity equation explicitly 
shows that mass diffusion is triggered by an unbalance of the local chemical potentials, so that equilibrium is attained 
whenever pA ^ Pb- So much for the continuum picture. 

For a finite value of r, an exact matching between momentum and continuity equations starting from continuum free- 
energy functional ([32]) is not so straightforward and more elaborate arguments are necessary. It is however possible 
to fix the bulk densities by introducing the following r-dependent functional 

(t) 

C-^Hpa, pb) = f^"\pA, pb) - 4^VpA ■ VpB - Xapa - XbPb (39) 

fi''^ {PA, Pb) = cIpa log PA + c%pb log Pb + c]^gABPAPB (40) 

where — |^ is the effective coupling renormalized by lattice discreteness effects (note that this is exactly the 

inverse of the reference density p'^g^ introduced earlier on). We note that in the long-time limit r ^ 1/2 {9{t) 1), 
we have 

(t-) 

9ab 9AB 

thus reproducing the continuum value. This r dependence of the effective coupling reflects into an analogue dependence 
of the bulk densities. The bulk minimization with respect of pA and ps, along the same lines as for the continuum 
case, leads to the following bulk equations (the same procedure leading to ([M]). with dxxPA,B = 0) 



dpA 



Cs9abPb = >^a ..y. 
-1W^ + 49abPa = Xb- 



The symmetry under the interchange pa ^ pB imposes — Xb- By subtracting the second from the first equation 
in (I4ip we obtain again the relation (j30p . thus showing that the manifold of {pA, Pb) minimizing the bulk free-energy 
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FIG. 3: Equilibrium bulk densities in the multicomponent Shan-Chen model with simple density-density repulsion 
(equation ([2| with i^J — 0, — 0). In this case we have chosen p'f^ = 1.6 in (|29p by setting qab ~ 0.345, 
po — 1 and T — 1.116071. We have then simulated a Id interface between two components at varying the total 
averaged density (p) in different numerical simulations. The numerical results are successfully compared with the 
prediction coming from the minimization of the free energy (|42p and well approximated by equations (|43|) . In the 
right figures we also show the typical profiles of the bulk free energies arising in the theory behind the simulations. 
For simplicity the bulk free energy has been normalized to minus a unit value at the two minima. All results are 
given in lattice Boltzmann units (LBU). 



is the same as the one obtained by imposing a zero diffusion current. Here again, in order to single out a point of 
minimum, we have to specify the total mass in the system, as stated in (j30p . The procedure gains transparency by 
replacing the densities with their their local sum (p — pA + Pb) and difference {(f) — pA — Pb)- The bulk- free energy 
functional takes then the following form: 

fi'\p. <!>) - iiP + ^) log ( + #(P - 0) log ( - #5i1(p^ - 0^) (42) 



It can be checked that, at a given value of p, this expression presents a double well structure, as soon as p > pg'^\ The 
two minima correspond to the two symmetric values of (j) (±0o) attained at equilibrium in the bulk phases. To be 
noted that the presence of the two minima reminds of the 'double tangent description' characterizing the minimization 
of a free-energy functional. In this simple case, due to the symmetric structure of the problem, we are left with a 
symmetric free energy and therefore the bulk densities can be directly extracted from those two minima. By Taylor 
expanding the full set of equations ([50)1 we obtain an analytical estimate of the solution for the two bulk densities 
(Pa B' Pa b where I, h stands for low and high density), namely: 



Pa.b = 



.Pa.b 



I ^-^=^^ ■ (43) 

(p) - y 30(pi,"' -6pi,"' V 45(p<"' )^ -10pi,"> (p) 



(r) 

This approach has been validated against numerical simulations. The results, referring to the case pg — 1.6, 
9AB — 0.345, Po = 1.0 and t = 1.116071 in lattice Boltzmann units (LBU), are shown in figured We have simulated 
a Id interface between two components at varying the total averaged density (p). The numerical results compare 
satisfactorily with the theoretical predictions based on the minimization of the free energy (|42p . In the right panel of 
the same figure, also shown are typical profiles of the bulk free energies arising in the numerical study. 



C. Multicomponent Model with Self-Interactions 



Having covered the case with purely repulsive inter-species interactions, we next consider the more general situation 
in which intra-species (self) interactions are included (equation Q with all interactions on). In this general case, the 
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condition of no mass diffusion ( Jq = 0) delivers: 



DAs{PA,PB)dxPs ^ ^ DBs{pA,PB)dxPs ^0 (44) 
s=A,B s=A,B 



with 



Daa = c%T (^^0{t) + ^{GaiPb-^a'^'a - 9abPbPa)^ (45) 

Dbb - cIt (^y 0(t) + ^(GsipA^B^'s - 9abPaPb)^ (46) 

with the usual symmetries: Dab = ^Dbb, Dba — —Daa- Also, a constant (Pg) pressure tensor across the interface 
is required: 

Pxx = {cIpa + cIpb + ^4^^!*^ + ^4^31*1 + clgABPAPB + 4uj - Txx + = Po (47) 

n = J:s=A,bGs2 ( 7(<9x*.)^ + I'i'sdxx'^s] + ^ (pAdxxPB + PbOxxPA + dxPAdxpB) (48) 



^2 

Txx - Y {G2Adx'fAdx'i'A + G2Bdx'i'Bdx'i'B + 9AB{dxPAdxPB + Q^PbOxPa)) (49) 

with the various effective couplings already defined in p?|) and Ki'x' defined in p3)) . These two 'conserved' currents 
must be matched with the total mass in the system. In the most general case, we expect two characteristic values of 
the sum of the two densities in the two bulks, corresponding to the four unknowns ^^41 Pa P^B' Pb- One can resort 
again to a minimization procedure based on the following free-energy density 

C{pa,Pb) - Mpa,Pb) + 4^|V*Ap + 4%^|V*s|' - 4^VpA ■ VpB - \aPa - XbPb (50) 
with the bulk contribution written as 

Gai f'^ ((,) Gbi [^'^ iO 

fb{pA,PB) =clpA\0gpA + clpB\0gpB + clgABPAPB + Cl^-PA J ^2 ^^ + 4"^^^ J ^2 

Note that, like in the purely repulsive case, this matches the equilibrium properties of our system in the limit t ^ 1/2, 
where lattice time discreteness can be ignored. Moreover, due to the presence of the pseudo-potentials '5, in order 
to make the Shan-Chen model compliant with such a kind of free energy, an extra-gradient term has to be added, 
as described in a recent paper Such extra-term is connected with variations of the pseudo-potentials across the 
interface and, at least for the case of a single-component fluid, it can be shown to be negligible to practical purposes. 
It is also worth noting that in the symmetric case Gai — Gbi (the one analyzed later in the paper) with the same 
pseudo-potential for both components ^ a = ^ b, we can use similar arguments as described in the previous subsection. 
In particular, we define the following r-dependent bulk free energy 

fi^\pA.PB)=4pA\0gpA + clpB\0gpB+C%g^llpAPB+4^PA ^^d^ + 4^PB (51) 

with G^]' = ^7y, ^^42 = ^-nd look for its (symmetric) minima. In figure [H we show the comparison between the 
results of minimization of this free energy and those by direct numerical simulations with the usual pseudo-potential 
^I^ — po(l — e^P/P°). The main parameters are po = 1-0, {p) = 1.75, gAB = 0.345 and t = 0.69 LBU, and different 
values of the self coupling parameters Gai = Ga2- Overall, satisfactory agreement is observed. 
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FIG. 4: Bulk densities in the numerical simulations with repusion plus pseudo-potentials (i.e. equation ([2| 
with _FJ — 0). Both pseudo-potentials are chosen in the same way, ^a.s ~ po(l — ^^p^.b/po-j fjxing po — 1.0, 
(p) = 1.75, r = 1.0, qab ~ 0.5785. By changing the self-coupling parameters in such a way that Gai ~ Ga2 we 
have computed the equilibrium bulk densities and also compared with the results coming from the minimization 
of the bulk (symmetric) free energy (|51|l . In the inset figure we also show the typical profile of the bulk free 
energy arising in the theory behind the simulations at Gia ~ —0.63. For simplicity the bulk free energy has been 
normalized to minus a unit value at the two minima. All results are reported in LBU. 



V. TRANSPORT PROPERTIES 



In this section we focus on the theoretical prediction of the surface tension of the two-component model. As 
previously discussed, this requires the correct identification of the off-diagonal component of the momentum-flux 
tensor Pab ■ For the sake of concreteness, we shall consider the simplest case of a one-dimensional stationary interface 
between the two fluids A and B. In view of equation (I14|) . the surface tension is given by 



CTAB = T^r^dx (52) 

J flat 

where /^j^jj is a short hand notation for integration across a flat interface separating the two fluids and developing 
across x. However, as pointed out by Shan & Chen |43i |, the time discretization induces an extra term on the r.h.s. of 
([5^ and, given its importance for the actual computation of the surface tension, in the following we shall generalize 
their treatment to the case of a two-component fluid. We start by writing the lattice kinetic equation for the total 
distribution function gi = fiA + fiB ■ 

gt{f+Ct,t+ 1) ~ g.,{f,t) = --[gt{f,t) - g[^''\pA, Pb,u,Fa,Fb)] 

T 

where the total equilibrium is simply the sum of the two single-component equilibria 

g'f'^ = /il"' {PA. u + tFa/pa) + fls'^ {Ps.u + tFb/pb). 
By unrolling the full expressions of f^'^^ and flj^\ we obtain: 



"S ^^S \ PA PB ^ 

where = -Fao + Fsa is the a-th component of the total force and 

vS^i^ = u + tFIp (54) 

is the total fluid velocity, including the shift due to the total force. To be noted that the term ( ''^'^p^^'' + "^""p^f^'' ~ 
-^^^) is missing in the original paper by Shan & Chen [i^, because these authors deal with a single-species fluid. 
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Next, following [4^, we estimate u'^'^'^^ by general considerations holding at steady state. For stationary solutions, we 
can assume no net mass transfer along any link connecting two lattice sites, which implies gi{r -\- Ci) = Q i(r ), where j 
is the mirror partner defined by the condition Cj — —Ci. Under this constraint, one derives the relation [43|: 

^2pu= --{pu~ pu^^'^y) (55) 
r 

which, combined with ([M)) . delivers: 

pu^^'"'' = j^r - F = r6'(r)F. 
This expression can then be used to evaluate the kinetic component of the pressure tensor pI7|) . 

i 

By assuming gi « g,-^''', the term giCiaCib delivers the following contribution: 

gi CiaCib = Cspdab +T (t) h T 1 (56) 

P \ PA PB P J 

where, as anticipated in the previous sections, we recognize the ideal gas equation of state, plus extra r-dependent 
contributions stemming from the forcing terms. This shows that discrete effects (both in time and space) introduce 
a correction to the surface tension, which must be taken into account in order to compute the value of aAB ■ We can 
now make use of the identity 

\ PA PB P ) p \ PA PB J \ PA PB ) ' 

Also, the condition of no mass-diffusion current p^ . gives: 

FAa FBa\ ^2nl^-^f ^aPA ^aPB A ^gg^ 



PA PB J \ PA PB 

Inserting (l57|) together with (l58|) into the rhs of (I56|l . finally delivers 

gt'^^CiaCih = clpSab + t'^0'^{t) 



FaFb ^ 4 PaPB ( daPA daPB \ ( dbPA 

— + "^s— 



PA PB J \ PA PB 



(59) 



which is precisely the result reported in ([Te 

In conclusion, the expression for the overall surface tension must take into account the contribution of both the 
potential energy and the (r-dependent) kinetic energy components of the pressure tensor: 



CTAB = 



T^,^dx+ / KilUx (60) 

flat J flat 



with Txx and K^'^J stemming from the interaction pressure tensor (equation (jl6[) ) and the r-dependent part of the 
kinetic pressure tensor (equation (jl3p ). respectively. 

The presence of the extra r dependent terms has been checked against numerical simulations with pure repulsion 
(equation ^ with FJ = 0, = ), as shown in figure [5] We have fixed p'f ^ ~ 0.91 in LBU and varied r in 
the simulations. The numerical results in figure [5] show the bare surface tension computed with and without the r 
corrections given in (|56p . as well as through the usual Laplace test, i.e.. by evaluating the difference between inner Pin 
and outer Pout equilibrium bulk pressure of two-dimensional droplets of radius i?, and extracting the surface tension 
from the Laplace's relation: 

P _ P - ^ 

^ in ^ out — • 

H 

The results clearly indicate that the correction terms are essential to achieve quantitative agreement with the Laplace's 
values. To be noted that the r-dependence of the equilibrium component of the kinetic pressure tensor, rhs of equation 
(I59p . which stems from the shifted velocity in the local equilibrium, disappears in the limit r — > 1/2. 
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FIG. 5: Surface tension in a model with pure particle-particle repulsion (equation ([2} with _FJ = 0, F° — 0) . We 
have fixed pj^' = 0.91 in LBU as defined in equation (|29|l and then varied r in the lattice Boltzmann simulations 
(all the numerical values reported are in LBU). Results show the bare surface tension computed with the simple 
interaction pressure tensor and also with the r corrections as described in (|56p . 



A. Achieving vanishingly-low surface tension for finite relaxation times 



Going back to the general expression of the forcing terms ^ it is interesting to observe that, once the values of the 
G-couplings in the full model are fixed, we can still tune the surface tension by suitably changing po in the model. 
For a fixed relaxation time r (say r = 1 LBU) this turns out to be a practical computational strategy to access the 
vanishing low surface-tension regime of interest for the simulation of micro-emulsions. Using the theory developed so 
far, we can now estimate the surface tension aAB as a function of the free parameter po appearing in equation 
Collecting the different terms coming from ((60|) . we obtain the following 



C^AB 



4 



dx 



flat 



Ga2 



G 



32 



Gai 

pI 



■d^PAd,pB+T^0'{T) 



PAPB ( dxPA dxPB 



PA 



PB 



(61) 



The exact computation of the integral in equation (j6ip requires the knowledge of the functions pa{x) and pb{x). 
However, useful insight can be gained by assuming that the sum of the two densities, pA + Pb = (p) is constant and 
that the leading contribution to the integral comes from the interface region, where pA ~ Pb- We can then expand 
about the point (/) = pa ^ Pb = that we consider located at the central point Xc- With these assumptions, we write 



PA~PB = 



2 



dxPA = -dxPB ~ dx(j)\x,, 



d^^A ~ e-^p'>'^P" O^paU dx^B « e'^p'>'^P" O^pbIx^. 
In this way, for r = 1 (LBU), equation (pTjl finally delivers 



(JAB ~ -Y'^AB{po){dy4>\xJ'^ 
with TiABipo) depending on the couplings and the parameter po as follows: 



■^AB 



iPo) 



{ 2Gab 2 
I Pi ip) 



- (Ga2 + GB2)e 



-(p>/Po 



(62) 



(63) 



and where 5w is the characteristic thickness of the interface. Equation (j63p shows that by increasing po the surface 
tension can be made negative, so that the condition Yiab{po) = stipulates a vanishing surface tension. Indeed, 
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FIG. 6: The function E(po) as defined in 163} for the set of parameters (p) = 1.23, Gab = 0.405, Ga2 = 9.1714, Gb2 = 8.45714. 
The crossover lies at about po = 0.72, as also evidenced in the inset. The line at zero surface tension is reported as a visual 
guideline. All results are reported in LBU. 
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FIG. 7: Surface tension as a function of po, as obtained by numerical simulations with the parameters (p) = 1.23, G'a = —15.0, 
Ga = 14.1, Gb = -14.0, Gb = 13.1, Gab = 0.405. This produces in Ga2 = 9.1714, Gb2 = 8.45714 in The crossover 

lies around po ~ 0.71 LBU, in close agreement with the analytical estimate. The line at zero surface tension is reported for 
eye-guiding purposes. All results are given in LBU. 



upon increasing p^, the positive contribution of repulsive interactions is weakened, whereas the negative contribution 
of self-interactions is enhanced, provided that Ga2 and Gb2 are both positive. In figure [6] we show the analytical 
computation of S^s as a function of po for the set of parameters (p) = 1.23, G\ = —15, = 14.1, = —14, 
G^ ~ 13.1, Gab = 0.405, corresponding to Ga2 — 9.17 and Gb2 = 8.46, all in LBU. The theory predicts a crossover 
of the surface tension to negative values at po ~ 0.72, quite close to the numerically observed result po ^ 0.717 (see 
figure [7]). This shows that the interplay between inter-species repulsion and intra-species repulsion/attraction is key 
to attain vanishing small values of the surface tension, which are in turn crucial to reproduce the physical properties 
described in the second part of this paper. 
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FIG. 8: Contours of the density field pA{x,y) obtained with a numerical simulation on a Nx x Ny — 512 x 512 grid. The 
parameters are those of the standard set at po — 0.83, as given in equations lf65|l . 



VI. NUMERICAL RESULTS 



Having discussed the major theoretical aspects of this model, we next proceed to present the results of numerical 
simulations. The baseline simulations are performed on a 2 dimensional grid Nx x Ny = 128 x 128, with occasional 
enlargements to Nx x Ny — 256 x 256 and N^ x Ny = 512 x 512. The two fluids are initialized with zero speed and 
random initial conditions for the two densities pA and pB- More specifically, we choose (pa) = (pb) = 0.612, with a 
standard deviation ±0.01 from the background density value. The couplings have been set to the following values in 
LBU: 



CX = -15.0, 
G% = -14.0, 
Gab = 0.405 



Gr 

A — 

Gr 

b — 



14.1 
13.1 



(64) 



defined as standard set at po = 0.7 and 



G'X = -9.0, 
G% = -8.0, 
Gab = 0.405 



Gr 



8.1 
7.1 



(65) 



defined as standard set at po — 0.83. The relaxation time is fixed to t = 1 (LBU), corresponding to a kinematic 
viscosity 1^ — 1/6 (LBU). The corresponding value of the surface tension is approximately aAB ^ 0.01 in both standard 
sets. The main difference between the two sets of parameter is that the standard set at po = 0.83 displays a more 
refined (in terms of computational grid points) interface. Moreover, the standard sets of parameters have been chosen 
in such a way that both components A and B are in the dense (liquid) phase. 



VII. FREE DYNAMICS OF THE DENSITY CONFIGURATION 



We begin by investigating the free configurational dynamics of the density field under the sole effect of internal 
interactions (no- forcing). The first observation is that, even after a very long time-span (hundreds of thousands 
time-steps) the fluid densities pAix, y) and pb{x, y) do not exhibit any macroscopic separation between the two fluids 
A and B. Instead, a multitude of metastable domains ("droplets") of fluid A in fluid B and viceversa is observed, 
as a result of the complex interplay between repulsive (short-range intcr-species and mid-range intra-species) and 
attractive (short-range intra-species) interactions. This is in line with other studies in solid state physics and soft 
matter [H, [63, [H, [63| . The final result is a rich configurational structure of the density field, as shown in figure [H 
The most salient feature of the density configurations is the formation of 'belts' of fluid A (B), entrapping bubbles of 
both fluids B and A inside. As we shall see shortly, these belts exert a major influence on the rheology of the fluid, 
and in particular, their formation/rupture is responsible for a number of features, such as dynamical heterogeneity 
and arrest, long-time relaxation, ageing effects and intermittency. 

The occurrence of belts of fluid A (B) entrapping fluid B (A), is well visible in figure [U where also shown (bottom 
panel) are the density cuts of species A, across the midline y/Ny = 0.5 for the two different standard sets of parameters 
at Po = 0.7 (see equations set (|64l) ) and po = 0.83 (see equations set (1551) ') . Although the details of the density contours 
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FIG. 9: Contours (top) and centerline cuts (bottom) of the density of the specie A, for the standard set of parameters at 
po ~ 0.7 (see equations set (|64pl and at po = 0.83 (see equations set (|65}). All results are reported in LBU. 



and profiles are clearly different in the two cases, the main qualitative feature, namely the presence of a multitude of 
metastable "droplets" of both fluids A and B, is well visible in both cases. Therefore, these "droplets" are naturally 
interpreted as the metastable structures which permit the two-fluid system to escape the fully-separated minimum- 
interface configuration. 



VIII. DYNAMIC RESPONSE UNDER APPLIED SHEAR 



In view of the rich morphology of the density field discussed in the previous section, it is natural to inspect the 
behavior of the two-fluid system under the effect of an external drive. To this purpose, we analyze the dynamic 
response to an externally applied shear flow of the form Ux{x,y) = Uosm{ky), Uy = 0, with fc = 1. This is realized 
by imposing a volumetric body force in the LB equation. The rheological properties of the fluid are measured by 
monitoring the following response function: 

R(t) = %M , ^ (66) 
Uo v{t) 

where U{k;t) is the Fourier transform of the line-averaged speed along the x direction, U{y;t) — U{x,y;t)/Nx, 
vq is the nominal kinematic viscosity of both fluids and P defines the effective viscosity of the two-fluid system. By 
construction, under undisturbed flow conditions, i? > 0, so that i? ^ 1 provides a direct measure of slowing-down 
through enhanced effective viscosity and eventually, structural arrest (R — 0). Baseline simulations are performed on 
a Nx X Ny ^ 128 x 128 grid, for up to 5 x 10^ LBU time steps. 



A. Cage formation and rupture 



A typical response function is shown in figure fTOl (lower panel), together with two snapshots of the density contours 
a.t t — 10^ and t — 3 10^ (upper panel). In the same figure, also shown is an indicator of the interface area (length in 
2d) between the two fluids, defined as follows: 

lAsit) = - ^V pA{x,y-t) - V pB{x,y;t). (67) 
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FIG. 10: Response function R{t) (squares) given in equation (|66|l and surface indicator lAsit) (circles) given in equation (|67p 
for the standard set of parameters at po = 0.7 (see equations set (|64|) '). The upper panels show two snapshots of the density 
field in a blocked and flowing state, respectively. Note that the flowing state is nonetheless characterized by a small fraction (a 
few percent) of the undisturbed flow speed, corresponding to values or R larger than 0. The dotted line at zero is reported for 
visual guidance. All results are given in LBU. 
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FIG. 11: The response function R{t) at three-different instants, t — 10^, 10"*, 5 10^ LBU, as a function of the applied forcing 
Uq. Simulations are carried out for the standard set of parameters at po = 0.7 (see equations set (|64p '). The inset reports the 
effective viscosity, i.e. the system average velocity versus the forcing amplitude. All results are reported in LBU. 



This figure provides a neat example of dynamical arrest (between t = h 10'' LBU and t ^ 1.5 10^ LBU, followed by 
a progressive recovery of the flow (from t ^ 1.5 10^ LBU to i ^ 3 10^ LBU, until the system starts to flow again, 
although with a 25-fold higher viscosity than the nominal one, i.e. R ^ 0.04 versus R—1. 

The two snapshots refer to a blocked configuration [t — 10^ LBU) and to a flowing one [t — i 10^ LBU), 
respectively. In the former, belts caging one fluid into another are well visible, which subsequently break down and 
disappear, thereby allowing the system to flow again. Consistently with this picture of cage rupture and annihilation, 
the interface length, as measured by Jab, is seen to decrease in going from the arrested to the cage-free flowing 
configuration. This picture clearly illustrates the vital role played by the cage structures on the global rheology of 
the two-fluid system. It is worth emphasizing that, due to the mesoscopic nature of the present model, the rupture 
of a single cage, corresponds to a large collection of atomistic events, and consequently it leads to observable effects 
on the overall rheology of the system. 

Next, we investigate the time dependence of the response function for different values of the shear forcing Uq. 

In figure nn we show a typical example for the response function R{t) at three-different instants, t = 10"^, lO'*, 5 10^ 
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LBU, and for different values of the forcing Uq. At short times the response is hnear with Uq for all investigated 
values of Uq, (Newtonian behavior). At longer times, however, a typical yield-stress threshold appears, i.e. the fluid 
starts to flow only beyond a critical value of the forcing, Uq ~ 0.03 LBU. Above this threshold, the fluid starts to 
flow at a higher rate (see also inset, reporting the effective viscosity) as compared to the short-time response, thereby 
providing evidence of non-newtonian, shear-thinning, behavior. 



B. Dynamics of correlations: ageing effects 




FIG. 12: Ageing of the system. Correlation function as defined in (|68p computed for different waiting times tw {tw ~ 5 10'', 
red squares, = 2 10^, green circles and = 3 10^, blue triangles ) with shear stress Uo = 0.02 LBU. The waiting time At„ 
is reported on a log scale. Simulations are carried out for the standard set of parameters at po — 0.7 (see equations set (|64|l '). 
The inset reports the correlation function for = 3 10^ and Uq = 0.03: with increasing shear stress the structural arrest 
disappears, as witnessed by a vanishing value of the correlation function in the limit At^ oo. All results are reported in 
LBU. 




FIG. 13: Ageing of the system at increasing shear stress. This configuration of the system represents a threshold-yield fluid. 
The system does not flow and shows ageing until a certain threshold of the applied forcing, above which the system decorrelates 
completely. In the figure 5* is the externally imposed shear and the waiting time Atw is reported on a log scale. All results are 
reported in LBU and the parameters are those of the standard set at po = 0.7 (see equations set (|64p '). 

We next inspect another typical phenomenon of soft-glassy matter, namely ageing. To this purpose, following upon 
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FIG. 14: Time evolution of the time-derivative of the average fluid velocity U{t) (inset), showing a characteristic intermittent 
behavior. In the main picture, shown is the probability distribution of the time-lapse (Ate) between two subsequent bursts 
of dU/dt. A typically power law decay ~ At~", with exponent a = 2 is observed. The figures refer to the standard set of 
parameters with po = 0.7 as given in equations set (|64|) and results are all expressed in LBU. 



the spin-glass literature [70| , we define the order parameter (f> = pA — Pb and compute its overlap, defined through 
the autocorrelation function: 

, <^(a;,y;t^)0(a;,y;t^ -f Ai^)) 
C{t^, At^) = — — (68) 

where is the waiting time, At^ is the time lapse between the two density configurations and brackets stand for 
averaging over an ensemble of reahzations. In figure 1121 we show the correlation function corresponding to three 
different waiting times, {t^ = 5 10'' LBU , red squares, = 2 10^ LBU, green circles and t^, = 3 10^ LBU, blue 
triangles), for a forcing amphtude Uq = 0.02 LBU. Ageing effects are clearly visible, in the form of a dependence 
of the time-decay of the correlation function on the waiting time i^,, and, more specifically, with an increasingly 
slower decay as the waiting time is increased. Moreover, the correlation function saturates to a non-zero value in 
the long-time limit (broken ergodicity), which is another typical signature of structural arrest (the system does not 
succeed to fully decorrelate) . This behavior shows qualitative changes upon increasing the forcing term. In the inset 
of the same figure, we show the correlation function for i„ = 3 10^ LBU and a slightly larger forcing, Uq = 0.03 LBU. 
With increasing shear stress, cages are broken, and the structural arrest disappears, thereby allowing the correlation 
function to decay to zero (see figure \T^ . The disappearance of structural arrest under sufficiently strong shear is 
again a distinctive feature of flowing soft-glassy materials (7l| and these results are in qualitative agreement with 
molecular dynamics simulations [t^ . 



C. Intermittency and Barkhausen noise 



Barkhausen noise is a well-known phenomenon displayed by disordered ferromagnetic samples under the effect of a 
slowly-changing magnetic field [73] . A small ramp- up in the magnetic field triggers one domain and the perturbation 
spreads to neighboring domains, producing an avalanche which results in a series of jumps in the magnetization, as 
the systems transits from one metastable state to another. Several experiments show that the distribution of size, 
duration and energy of the Barkhausen jumps exhibit a power-law decay. The present two-fluid model also shows 
evidence of Barkhausen- like intermittency in the time-derivative of the response function. In figure 1141 we show the 
probability distribution of the time-lapse Ate between subsequent bursts (also called 'events') of the response function 
(see inset). Interestingly, such distribution follows a power-law distribution ~ At~", with a ~ 2. This invites a further 
analogy between the fluid cages discussed previously and the magnetic domains responsible for Barkhausen effects 
in disordered ferromagnets. The systematic exploration of the dependence of these Barkhausen-like effects on the 
various parameters of our system, is left as an interesting topic for future research. 
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IX. SENSITIVITY TO INITIAL CONDITIONS 



In the previous sections we have iUustrated several phenomena typical of soft-glassy materials. A natural question 
arises as to the degree of robustness of these phenomena towards changes in the initial conditions and size of the 
system. Although a systematic exploration of these effects requires a study of its own, in the following we provide 
some preliminary information. As expected, the detailed dynamics of the response function shows a strong sensitivity 
to the noise realization, with some configurations reaching a plateau in the early stage of the evolution (see figure 
I15|l . while others never attaining any plateau within the entire simulation span. In order to probe the robustness of 
the response function R{t) towards changes in the random realization of the initial conditions, we have performed a 
series of 100 simulations by changing the noise realization at a fixed variance of the initial density. Notwithstanding 
the qualitative differences in the detailed response function, the main picture portrayed in the previous sections, 
namely arrested flow due to formation of fluid cages, and restored flow upon cage rupture, is found to apply to all 
simulations. To better appreciate the statistical dynamics of the present system, in figure \W[ top panel, we show the 
time evolution of the Kurtosis IC{t) = ^^;^(^/)2^2 of the response R{t), as computed from the set of 100 realizations. This 
figure shows clear evidence of large fluctuations in the first half of the evolution, followed by a more quiescent stage 
in the second half. To be noted that, even in the quiescent stage, the Kurtosis is still around IC{t) ^ 5, hence well 
above the Gaussian value JC(t) = 3, thereby confirming the strongly fluctuating nature of the phenomenon. A similar 
message is conveyed by the bottom panel of the same figure, which reports the average value {R(t)), along with the 
variance, as a function of time. From this figure, we see that the variance is generally comparable to the mean value, 
sometimes even larger. The intermittent nature of the response R{t) is further highlighted in figure fT7| which shows 
the probability distribution function of R{t), sampled over three close-by time-slices. This pdf exhibits intermittent 
tails on both negative and positive sides, with a slight prevalence of the latter, consistently with the positive sign of 
(i?(t)). 



X. SUMMARY AND OUTLOOK 



Summarizing, we have provided a theoretical analysis of a two-component lattice Boltzmann model with mid- 
range intra-molecular repulsion and short-range inter-molecular repulsion. In particular, equilibrium densities and 
the surface tension as a function of the main parameters of the model, have been computed and shown to exhibit 
satisfactory agreement with numerical tests. We have also presented a series of numerical simulations proving the 
capability of this system of reproducing many distinctive features of soft material behavior, such as slow-relaxation, 
anomalous enhanced viscosity, caging effects, aging under shear and Barkhausen intermittency. The present lattice 
kinetic model caters for this very rich physical picture at a computational cost only marginally exceeding the one for 
a simple fluid. As a result, it should be possible to use it for future investigations of the non-equilibrium rheology. 
In particular, it may be useful to get new insights in the coexistence of liquid and solid regions (shear localization, 
shear banding, cracks) as observed with emulsions [tiI, [tII , foams [tI, [tI, [Tg , worm- like micelles J8„ 79] and granular 
materials [sol [8l|. Still, such a hydro-kinetic method might be interesting to treat the issue of dilatancy in foams 
observed in recent experiments [84] • In order to analyze those systems, on going research is devoted to a systematic 
investigation of the system behavior at different concentrations of the two species, its sensitivity to initial conditions 
and flnite-size effects, as well as its response to time-dependent loads. Also of current interest are extensions to 
three-component fluids, in order to account for the explicit presence of surfactants [s^. 



Appendix: heuristic mapping to physical units 



One of main advantages of the present mesoscopic approach is to provide access to hydrodynamic scales at an 
affordable computational cost. In order to appreciate this point, it is of interest to discuss the conversion between 
LB and physical units. The spatial units, namely the LB spacing Ax, can be estimated by fixing the surface tension 
according to the following relation (subscript phys denotes physical units): 

kT 

CTphys ~ CTLBJ^^ (69) 

where the subscript LB denotes the value in LBU. For micro-emulsions, we may estimate aphys ^ 10^* N/m, so that 
at standard conditions (T = 300°), a LB surface tension aLB ~ 0.01 corresponds to Ax = . /kT '^^'^ ~ 10~^ m. 
This means that a x Ny — 128 x 128 simulation covers a squarelet of about 0.1 micron in side. Similarly, the time 
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FIG. 15; Response function R{t) (squares) given in equation (|66|l and surface indicator lAsit) (circles) given in equation (|67p 
for two different realizations of the initial conditions, obtained by changing the seed of the random number generator. The 
runs are performed with the standard set of parameters at po = 0.7 as given in equations set (|64|l . The dotted line at zero is 
reported as a visual guidance. All results are reported in LBU. 



units (the LB time step At) can be estimated by fixing the kinematic viscosity according to the relation: 

(Ax)^ 

Vphys « VLB (70) 

By taking VpUys ^ 10"^ m? / s and vlb ^ 0.1, a lattice spacing Ax ~ 10~^ to, would yield /Si ^ 10"^'^ s. As a result, 
a 10^ time-step simulation covers about 0.1 ^s. These values are only marginally higher than those typically used in 
Molecular Dynamics simulations. However, the point is that the present model lends itself to substantial upscaling 
both in space and time, while still presenting an affordable computational cost. For instance, preliminary simulations 
on a X Ny = 1024 x 1024 grid, span 10^ lattice time steps in about one-day elapsed time on Graphical Processing 
Units architecture j83 |. Such simulations cover a square domain about some microns in side, over a time span of some 
microseconds, a way beyond the capabilities of standard Molecular Dynamics or Monte Carlo simulations. 
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FIG. 16: Evaluation of the response function R{t) under the apphcation of a shear forcing Uo kept fixed in the numerical 
simulations. We have chosen 200 equispaced times between t = and t = 10® (LBU) and run 100 numerical simulations by 

changing the random initial conditions. The upper plot shows the Kurtosis, IC{t) = j^^^j^z, where (...) refers to the average 
over the various numerical runs at a fixed t. Intermittency is clearly visible from this plot. For the sake of clarity, we have 
reported the value of the Kurtosis for a Gaussian variable (equal to 3, dotted line). In the lower plot, the response function and 
the variance (errorbars) is then evaluated and plotted as a function of time. The line at zero is reported as a visual guideline. 
For the numerical simulations, we have used the standard set of parameters at po = 0.7 as given in equations set (|64p . with a 
shear forcing term Uo/N^ = 0.1/128 in LBU. The resolution is x Ny ^ 128 x 128. 
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FIG. 17: Plot of the probability density function for the response function R{t) under the application of a shear forcing Uo/Nx, 
kept fixed in the numerical simulations. The probability has been obtained from the analysis of the response function at 
three distinct but very close times t = 4.8 10^,4.9 10^,5 10^ (LBU) in 100 numerical simulations with different random initial 
conditions (the choice of 3 close instants is meant to enhance the statistics) . The numerical simulations are performed with the 
standard set of parameters at po — 0.7, as given in equations set (|64|) with a shear forcing term Uo/N^ ~ 0.1/128 in LBU. The 
resolution is x Ny = 128 x 128. 
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